clear all
tic
mcvar=zeros(51,1);
for ii=1:10
    randseed=randn(5000,5000);
    sampleW=oupath(randseed,0,1,1,0);
    %sampleOU=oupath(randseed,1,1,1,1);
    mcvar=mcvar+sum(midpoint(sampleW.^2).^2,2);
    %mcvar=mcvar+sum(sampleW2.*sampleOU2,2);
end
mcvar=mcvar/50000;
toc

% a=1;b=1;m=1;yic=1;
time=(0:50)*0.5/5000;
% exactvar=b^2/2/a*(1-exp(-2*a*time))'+(m+(yic-m)*exp(-a*time)).^2';
plot(time,mcvar)
% hold on
% plot(time,exactvar,'r')